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Abstract 

We generalize the recent invariant polytope algorithm for computing the joint spec¬ 
tral radius and extend it to a wider class of matrix sets. This, in particular, makes the 
algorithm applicable to sets of matrices that have finitely many spectrum maximizing 
products. A criterion of convergence of the algorithm is proved. 

As an application we solve two challenging computational open problems. First we 
find the regularity of the Butterfly subdivision scheme for various parameters u. In the 
“most regular” case w = ^, we prove that the limit function has Holder exponent 2 
and its derivative is “almost Lipschitz” with logarithmic factor 2. Second we compute 
the Holder exponent of Daubechies wavelets of high order. 

Keywords: matrix, joint spectral radius, invariant poly tope algorithm, dominant prod¬ 
ucts, balancing, subdivision schemes, Butterfly scheme, Daubechies wavelets 


1 Introduction 

The joint spectral radius of a set of matrices (or linear operators) originated in early sixties 
with Rota and Strang |3l] and found countless applications in functional analysis, dynamical 
systems, wavelets, combinatorics, number theory, automata, formal languages, etc. (see 
bibliography in [HI [T5, EJ (2Sj)- We focus on finite sets of matrices, although all the results 
are extended to arbitrary compact sets. If the converse is not stated, we assume a fixed 
basis in W 1 and identify an operator with the corresponding matrix. Everywhere below 
A = {Ri,..., A m } is an arbitrary family of d x d-matrices, A k is the set of all m k products 
of k matrices from A (with repetitions permitted). A product II G A k , n G N, is called 
simple if it is not a power of a shorter product. 
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Definition 1 The joint spectral radius of a family A is [3^] 

p(A) = lim max ||A|| 1//fc . (1) 

k-to o A&A k 

This limit exists and does not depend on the matrix norm. In case m — 1 (i.e., A = {Ai}), 
according to Gelfand’s theorem, the joint spectral radius is lini/ ; .^ 00 ||>1*|| 1 / fc , i.e., coincides 
with the usual spectral radius p{Af), which is the maximal modulus of eigenvalues of A\ 
(see, for instance PI EGO)- 

The joint spectral radius has the following geometrical meaning: p(A) is the infimum of 
numbers v for which there is a norm || • \\ v in W 1 such that in the induced operator norm, 
we have || A|| v < v, i — 1,..., m. In particular, p(A) < 1 if and only if all operators from A 
are contractions in some (common) norm. Thus, the joint spectral radius is the indicator of 
simultaneous contractivity of operators Ai, ..., A m . 

Another interpretation is due to the discrete dynamical system: 

x(k + 1) = A{k)x(k ), k — 0,1, 2,..., 

where each matrix A(k) is chosen from A independently for every k and x(0) G \ {0}. 
Then p(A) is the exponent of fastest possible growth of trajectories of the system: the 
maximal upper limit lim sup log among all the trajectories {xk}k >o is equal to logp(M). 

k —^OO 

In particular, p(A) < 1 if and only if the system is stable, i.e., x(k) —> 0 as k —> oo, for every 
trajectory. 

The problem of computing or estimating the joint spectral radius is notoriously hard. 
This is natural in view of the negative complexity results of Blondel and Tsitsiklis 13 d. 
Several methods of approximate computation were elaborated in the literature PI SI El El 
ITT1 Ifil iT5l 25] ESI EHs SZ]- They work well in low dimensions (mostly, not exceeding 5 — 8). 
When the dimension growth, then ether the estimation becomes rough or the the running 
time grows dramatically. In recent work na we derived the invariant polytope algorithm that 
finds the exact value of p(A) for the vast majority of matrix families in dimensions up to 20. 
For sets of nonnegative matrices it works faster and finds the exact value in dimensions up 
to 100 and higher. We conjectured in mi that the set of matrix families A for which this 
algorithm finds p(A) within finite time is of full Lebesgue measure. Several open problems 
from applications have been solved by using this method. However, it cannot handle one 
important case, which often emerges in practice: the case of several spectrum maximizing 
products. In this paper we generalize this method making it applicable for a wider class of 
matrix families, including that special case. To formulate the problem we need some more 
facts and notation. 

The following double inequality for the joint spectral radius p = p(A) is well known: 

max[p(A)] 1 / fc < p < maxllAH 1 ^, k } t e N. (2) 

AeA k AeA i 

Moreover, both parts of this inequality converge to p as k —> oo. In fact, 

lim max ||A|| 1/,fc = p (by definition) and lim sup max[p(A)] 1//fc = p (see m- All algorithms 
fc-s>oo AeA k k —^oo A&A k 
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of approximate computation of the joint spectral radius are based on this inequality. First, 
one finds the greatest value of max^g^t [p(A)] 1//fe (the left hand side of (J2D) among all reason¬ 
ably small k, then one minimizes the value max^k UAH 1 /* (the right hand side of (J2D) by 
choosing an appropriate matrix norm || • ||. Thus maximizing the lower bound and minimiz¬ 
ing the upper one we approximate the joint spectral radius. Sometimes those two bounds 
meet each other, which gives the exact value of p. This happens when one finds the spectrum 
maximizing product A = II G A k and the norm || • || for which both inequalities in ([2]) become 
equalities. 

Definition 2 A simple product II G A n is called the spectrum maximizing product (s.m.p.) 
if the value [p{Tl)] 1 / n is maximal among all products of matrices from A of all lengths n G N. 

Let us remark that an s.m.p. maximizes the value [p(II)] 1 / n among all products of our 
matrices, not just among products of length n. Observe that for any s.m.p., II G A n , we 
have [p(II)] 1 / ri = p(A). Indeed, from ([2]) it follows that [p(n)] 1//n < p(A). If this inequality 
is strict, then there are k G N such that max^g^fcf/^A)] 1 /^ > [p(n)] 1,/n (because of the 
convergence property), which contradicts to the maximality of II. Thus, to find the joint 
spectral radius it suffices to prove that a given product II G A n is an s.m.p.. The invariant 
polytope algorithm mi proves the s.m.p. property of a chosen product II by recursive 
construction of a polytope (or more general P C C d , although here we consider for simplicity 
real polytopes) P Ci d such that AiP C [p(II)] 1//n P, i = 1, ..., m. 

In the Minkowski norm || • ||p generated in by this polytope, we have max^g^ ||A||p < 
[p(II)] 1 / n and || - ||p is said to be an extremal norm for A. Applying ()2jl for k — n (left hand 
side inequality) and for t = 1 (right hand side) we conclude that [p(II)] 1 / n = p{A). 

Note that if a product II G A n is s.m.p., then so is each of n its cyclic permutations. If 
there are no other s.m.p., then we say that the s.m.p. is unique meaning that it is unique 
up to cyclic permutations. 

The disadvantage of the polytope algorithm is that it is guaranteed to work only if the 
family A has a unique s.m.p. Otherwise the algorithm may not be able to construct the 
desired polytope within finite time, even if this exists. The uniqueness of s.m.p. condition, 
although believed to be generic, is not satisfied in many practical cases. For example, it 
happens often that several matrices A % G A are s.m.p. (of length 1). The extension of our 
algorithm presented in this paper works with an arbitrary (finite) number of s.m.p.’s. We 
prove the theoretical criterion of convergence of the algorithm and apply it to solve two long 
standing open problems: computing the Holder regularity of the Butterfly subdivision scheme 
(Section O and computing the regularity of Daubechies wavelets of high order (Section [5]). 


2 Statement of the problem 

We begin with a short description of the invariant polytope algorithm from HI for computing 
the joint spectral radius and finding an extremal polytope norm of a given family A = 
{Ai,..., A m }. We make a usual assumption that the family is irreducible, i.e., its matrices 
do not have a common nontrivial invariant subspace. Recall that for a reducible family the 
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computation of the joint, spectral radius is obtained by solving several problems of smaller 
dimensions up. For a given set M C we denote by co(M) the convex hull of M and 
by absco(M) = co{M, — M} the symmetrized convex hull. The sign x denotes as usual the 
asymptopic equivalence of two values (i.e., equivalence up to multiplication by a constant). 

The invariant poly tope algorithm (see [30 1 [ISl HTj ). 

Initialization. First, we fix some number n and find a simple product II = A dn ... A dl 
with the maximal value [p(II)] 1 / n among all products of lengths n <n. We call this product 
a candidate s.m.p. and try to prove that it is actually an s.m.p. Denote p c = [p{Ii)} l ^ n and 
normalize all the matrices Ai as A * = p~ 1 A i . Thus we obtain the family A and the product 
II = A dn ... A dl such that p(II) = 1. For the sake of simplicity we assume that the largest 
by modulo eigenvalue of II is real, in which case it is ±1. We assume it is 1, the case of 
—1 is considered in the same way. The eigenvector corresponding to this eigenvalue 
is called leading eigenvector. The vectors = A dj _ 1 ■ ■ ■ A dl v^ l \ j = 2,..., n, are leading 
eigenvectors of cyclic permutations of IT. The set P = ... , is called root. Then 

we construct a sequence of finite sets V,; C and their subsets 7 Zi C V t as follows: 

Zero iteration. We set Vo = IZo = PL. 

k -th iteration, k > 1 . We have finite set Vk-i and its subset TZk-i- We set Vk = 
Vk- 1 , 7 Zk = 0 and for every v G TZk-i, A G A, check whether Av is an interior point 
of absco(VA.) (this is an LP problem). If so, we omit this point and take the next pair 
(v, *4.) G TZk-i x A, otherwise we add Av to Vk and to 7 Zk- When all pairs (v, A) are 
exhausted, both Vk and IZk are constructed. Let Pk = absco(W). We have 

Vk = Vk- 1 U 7 Zk , Pk = co {AiP k -i, ..., A m P k -i} . 

Termination. The algorithm halts when Vk = Vfc_i, i.e., TZk = 0 (no new vertices are 
added in the 7c-th iteration). In this case Pk-i = Pk, and hence Pk-i is an invariant polytope, 
II is an s.m.p., and p{A) = [p(II)] 1//n . End of the algorithm. 

Actually, the algorithm works with the sets Vk only, the polytopes Pk are needed to 
illustrate the idea. Thus, in each iteration of the algorithm, we construct a polytope Pk C 
store all its vertices in the set Vk and spot the set 7 Zk C Vk of newly appeared (after the 
previous iteration) vertices. Every time we check whether APk C Pk- If so, then Pk is an 
invariant polytope, ||Aj||p fc < 1 for all i, where || • ||p fc is the Miknowski norm associated to 
the polytope Pk, and II is an s.m.p. Otherwise, we update the sets Vk and 7 Zk and continue. 

If the algorithm terminates within finite time, then it proves that the chosen candidate 
is indeed an s.m.p. and gives the corresponding polytope norm. Although there are simple 
examples of matrix families for which the algorithm does not terminate, we believe that 
such cases are rare in practice. In fact, in all numerical experiments made with randomly 
generated matrices and with matrices from applications, the algorithm did terminate in finite 
time providing an invariant polytope. The only special case when it does not work is when 
there are several different s.m.p. (up to cyclic permutations). In this case the algorithm 
never converges as it follows from the criterion proved in mi- The criterion uses the notion 
of dominant product which is a strengthening of the s.m.p. property. A product II G A n is 
called dominant for the family A if there is a constant 7 < 1 such that the spectral radius 
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of each product of matrices from the normalized family A = [pilPf^^A, which is neither a 
power of fl nor one of its cyclic permutations, is smaller than 7 . A dominant product is an 
s.m.p., but, in general, not vice versa. 

Theorem A [TJ]. For a given set of matrices and for a given initial product n , the invariant 
polytope algorithm terminates within finite time if and only if II is dominant and its leading 
eigenvalue is unique and simple. 

Note that if there is another s.m.p., which is neither a power of fl nor of its cyclic 
permutation, then fl is not dominant. Therefore, from Theorem A we conclude 

Corollary 1 If a family A has more than one s.m.p., apart from taking powers or cyclic 
permutations, then, for every initial product, the invariant polytope algorithm does not ter¬ 
minate within finite time. 


The problem occurs in the situation when a family has several s.m.p., although not generic, 
but possible in some relevant applications. Mostly those are s.m.p. of length 1, i.e., some 
of matrices of the family A have the same spectral radius and dominate the others. This 
happens, for instance, for transition matrices of refinement equations and wavelets (see 
Sections [5l [ 6 ]). In the next section we show that the algorithm can be modified and extended 
to families with finitely many spectrum maximizing products. 

Let a family A have r candidate s.m.p.’s hi!,..., II,., r > 2. These products are assumed 
to be simple and different (up to cyclic permutations). Denote by 77 the length of flj. Thus, 
[p(n 1 )] 1 / ni = . . . = [p( II r )] 1 / nr ' = p c . Let A = pf l A be the normalized family, 77 be a leading 
eigenvector of LI* (it is assumed to be real) and "H* = • • •, } be the corresponding 

roots, i — 1..., r. The first idea is to start constructing the invariant polytope with all roots 
simultaneously, i.e., with the initial set of vertices 


Vo 


u un, 



1 , . . . , Hi , 7 = 1 , 



However, this algorithm may fail to converge as the following example demonstrates. 

Example 1 Let Ai and A 2 be operators in M 2 : Ai is a contraction with factor | towards 
the OX axis along the vector (1, 4) T , A 2 is a contraction with factor | towards the OY axis 
along the vector (1, — 2 ) 7 . The matrices of these operators are 



Clearly, both A 1 and A 2 have a unique simple leading eigenvalue 1; V\ = (1, 0) T is the leading 
eigenvector of Ai and v 2 = (0,1) T is the leading eigenvector of A 2 . 

The algorithm with two candidate s.m.p.’s Hi = Ai,n 2 = A 2 and with initial vertices 
Vo = {vi,v 2 } does not converge. Indeed, the set absco(Vfc) has an extreme point A%v 1 , which 
tends to the point 2v 2 as k —> 00 , but never reaches it. 

On the other hand, the same algorithm with initial vertices Vo = {ui, 3 t> 2 } terminates 
immediately after the first iteration with the invariant polytope (rhombus) P = absco(Vo) = 
co{±ui, ± 3 u 2 }- Indeed, one can easily check that AiP C P, i = 1,2. Therefore, A± and A 2 
are both s.m.p. and p(A) = 1 . 
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This example shows that if the algorithm does not converge with the leading eigenvectors 
Vi,. .., v r (or with the roots Tfi,..., Fir ), it, nevertheless, may converge if one multiplies these 
eigenvectors (or the roots) by some numbers aq,..., ay. In ExampledJwe have aq = 1 ,a 2 = 3. 
We call a vector of positive multipliers a = (op,..., a r ) balancing vector. 

Thus, if a family has several candidate s.nr.p.’s, then one can balance its leading eigenvec¬ 
tors (or its roots) i.e., multiply them by the entries {ai}( =1 of some balancing vector a > 0 
and start the invariant polytope algorithm. In the next section we prove that the balancing 
vector a for which the algorithm converges does exist and can be efficiently found, provided 
all the products Hi,..., IR are dominant (meaning the natural extension of dominance from 
a single product to a set of products). Thus, the corresponding extension of the invariant 
polytope algorithm is given by Algorithm [H 

A crucial point in Algorithm [T] is Step 5, where suitable scaling factors aq,..., ay have to 
be given. Then Algorithm [ 1 ] essentially repeats the invariant polytope algorithm replacing 
the root FL by the union of roots ot\FL\, ..., a r FL r - Finding the scaling factors that provide 
the convergence of the algorithm is a nontrivial problem. A proper methodology to compute 
them in an optimal way (in other words, to balance the leading eigenvectors) is derived in 
the next section. 

Remark 1 Till now we have always assumed that the leading eigenvalue is real. This is not 
a restriction, because the invariant polytope algorithm is generalized to the complex case as 
well (Algorithm C in [T71 U51 122] )- For the sake of simplicity, in this paper we consider only 
the case of real leading eigenvalue. 


3 Balancing leading eigenvectors. The main results 

3.1 Definitions, notation, and auxiliary facts 

Let ILi,..., II r be some products of matrices from A of lengths ni,... ,n r respectively. They 
are assumed to be simple and different up to cyclic permutations. We also assume that all 
those products are candidates s.nr.p.’s, in particular, [p(Ili)] 1//ni = ... = [p(n r )] 1 / nr = p c . 
We set A = pf x A and denote as M the supremum of norms of all products of matrices from 
A. Since A is irreducible and p(A) = 1, this supremum is finite [3]. By A* = {A\, ..., A^} 
we denote the family of adjoint matrices, the definition of A* is analogous. 

To each product A bn ... A bl we associate the word b n ... bi of the alphabet {1 ,..., m}. 
By ab we denote concatenation of words a and b, in particular, a n = a.. .a (n times); the 
length of the word a is denoted by |a|. A word is simple if it is not a power of a shorter 
word. In the sequel we identify words with corresponding products of matrices from A. 

Assumption 1 We now make the main assumption: each product Lb has a real leading 
eigenvalue (RemarkU)) , which is either 1 or —1 in this case. For the sake of simplicity, we 
assume that all A* = 1 (the case A* = —1 is considered in the same way). We denote by v t 
the corresponding leading eigenvector of II; (one of them, if it is not unique). 
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Algorithm 1: The invariant polytope algorithm extension 


1 

2 

3 

4 

5 

6 

7 

8 
9 

10 

11 

12 

13 

14 

15 

16 

17 

18 
19 


Data. A {A},..., A m j-, k max 

Result: The invariant polytope P , spectrum maximizing products, the joint spectral 
radius p(A) 

begin 

Compute a set of candidate spectrum maximizing products IR,..., II r ; 

Set p c := p( IR ) 1 /™ 1 and A := pj L A; 

Compute v\,..., v r , leading eigenvectors of IR,..., IR with ||n,-|| = 1 for all j; 
Form the roots TR,..., 7R; 

Provide the positive scaling factors ati,...,a r < 1; 

Set Vo := {njHjj'j , H 0 := V 0 ; 

Set k — 1; 

Set E = 0; 

while E = 0 and k < k max do 
Set V k = V fc _i, TZ k = 0; 
for u G Tlk-ij and for z = 1 ,..., m do 
if AiV G int(absco(Vfc)) then 
Leave Vfc,7 Z k as they are; 
else 

|_ Set V fc := V fc U {Ru}, 7^ fc := TZ k U {Ru}; 

if lZ k = 0 then 

Set E = 1 (the algorithm halts) ; 
else 

Set k k + 1 ; 

if E = 1 then 

return P := absco(Vfc) is an invariant polytope; 

III,... II r are s.m.p.; 

p(A) = p c is the joint spectral radius; 

else 

print Maximum number of iterations reached; 
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Clearly, the corresponding adjoint matrix IT* also has a leading eigenvalue 1, and a real 
leading eigenvector v*. If the the leading eigenvalue is unique and simple, then (v*,vf) ^ 0. 
In this case we normalize the adjoint leading eigenvector v* by the condition (v*,vf) = 1. 

Take arbitrary i — 1,..., r and consider the product 11* = Ad n ■ ■ ■ A^, where n = n* and 
d s = (for the sake of simplicity we omit the indices i). The vectors v\ 1 ^ = u*, = 

A^Vi, ... , vl' >] = Ad, nl • • ■ A^Vi are the leading eigenvectors of cyclic permutations of the 
product Ilj. The set Hi = {v* * \ ..., vf } } is a root of the tree from which the polytope 
algorithm starts. Let P*y = absco { A p H t | p = 0,..., k } be the polytope produced after 
the fc-tfi iteration of the algorithm started with the product 11*, or, which is the same, with 
the root Hi- This polytope is a symmetrized convex hull of the set V*y of all alive vertices 
of the tree on the first k levels. In particular, V* j0 = 'Hi and P* j0 = absco('H.j). We denote 
by Pi, oo and V* )0O the union of the corresponding sets for all k G Nil {0}. If Algorithm Q] 
terminates within finite time, then V* j00 is finite and P* !00 is a polytope. If 11* is an s.m.p., 
i.e, if p(LI.j) = 1, then, by the irreducibility, the set P ioo is bounded [3J. 

Now assume that all II* have unique simple leading eigenvalues, in which case all the 
adjoint leading eigenvectors v* are normalized by the condition ( v*,Vi ) = 1. For an arbitrary 
pair (i,j) 6 {1,..., r} 2 , and for arbitrary k G {0} U N U {oo}, we denote 

Qif = sup | (v*, z) |, (3) 

(k) 

Thus, q\J is the length of the orthogonal projection of the convex body P*y onto the vec¬ 
tor v*. In particular, qf) = max ze% |(n*, z)\ = max i=lr . i!n . |(n*, u®)|. Sometimes we omit 
the superscript k if its value is specified. For given j e {1,..., r} and for a point x G R d , we 
denote 

qf\x ) = _ max | (v *, ftx) |; q^x) = qf°\x) = sup | (v*, ftx) | (4) 

fieAp,p=o,...,k fi£Ap,p>o 

Note that for x = n*, the sets {Ida; | II G A p , p = 0,..., k} and V*y have the same sym¬ 
metrized convex hull P*^. Therefore, the maxima of the function | ( v *, z ) | over these two 
sets coincide with the maximum over P*^. Hence, comparing (|3]) and ([4|) gives 

Qj k \vi) = 4 fe) > k - ° : ^( v i) = 

For an arbitrary balancing vector a = (aq,..., ay), we write a -H = ..., a r H r }. Our 

aim is to hnd a balancing vector such that the polytope algorithm starting simultaneously 
with the roots aH terminates within finite time. 

Definition 3 Let k G {0} U N U {oo}. A balancing vector a G Ik} is called k-admissible, if 

a i Qij < a i h J = ■ ■ ■, r, if- j . (5) 


An oo-admissible vector is called admissible. 


(k) 

Since the value gh ; is non-decreasing in k, we see that the ^-admissibility for some k implies 
the same holds true for all smaller k. In particular, an admissible vector is fc-admissible for 
all k > 0 . 

We begin with two auxiliary facts needed in the proofs of our main results. 

Lemma 1 If a dx d matrix A and a vector x E are such that || Ax — x|| < e ||x||, 

then A has an eigenvalue A E C for which | A — 11 < C(d) ||A|| e 1 ' d , where C(d ) depends only 
on the dimension d. 

See pH>] for the proof. The following combinatorial fact is well-known: 

Lemma 2 Let a, h be two simple nonempty words of a finite alphabet, n, k > 2 be natural 
numbers. If the word a n contains a subword b k such that \b k \ > |a|, then b is a cyclic 
permutation of a. 

Now we extend the key property of dominance to a set of candidate s.m.p.’s. 

Definition 4 Products IIi,...,n r are called dominant for the family A if all numbers 
WH,)]■/”., a are equal (denote this value by p c ) and there is a constant 7 < 1 such 

that the spectral radius of each product of matrices from the normalized family A = pf l A 
which is neither a power of some hf nor that of its cyclic permutation is smaller than 7. 

3.2 Criterion of convergence of Algorithm [T] 

If the products Ili,...,n r are dominant, then they are s.m.p., but, in general, not vice 
versa. The s.m.p. property means that the function /(II) = [p(II)] 1//?l (77, is the length of II) 
defined on the set of products of the normalized family A attains its maximum (equal to one) 
at the products IT, and at their powers and cyclic permutations. The dominance property 
means, in addition, that for all other products, this function is smaller than some 7 < 1 . 
This property seems to be too strong, however, the following theorem shows that it is rather 
general. 

Theorem 1 Algorithm |7] with the initial products III,..., n r and with a balancing vector a 
terminates within finite time if and only if these products are all dominant, their leading 
eigenvalues are unique and simple and a is admissible. 

The proof is in Appendix 1. 

Remark 2 Theorem [T] implies that if the algorithm terminates within finite time, then the 
leading eigenvalues of products IT, must be unique and simple. That is why we defined 
admissible balancing vectors for this case only. 

If Algorithm Q] produces an invariant polytope, then our candidate s.m.p.’s are not only 
s.m.p.’s but also dominant products. A number of numerical experiments suggests that the 
situation when the algorithm terminates within finite time (and hence, there are dominant 
products) should be generic. 
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3.3 The existence of an admissible balancing vector 

By Theorem [TJ if all our candidate s.m.p.’s {ITj}[ =1 are dominant and have unique simple 
leading eigenvalues, then balancing the corresponding roots by weights oq,... ,a r 

we run the algorithm and construct an invariant polytope, provided the balancing vector 
a is admissible. A natural question arises if an admissible vector always exists. The next 
theorem gives an affirmative answer. Before we formulate it, we need an auxiliary result. 

Lemma 3 For given coefficients qff the system ([5]) has a solution a > 0 is and only if for 
every nontrivial cycle (A,... ,i n ) on the set {!,..., r}, we have (with i n+ 1 = i\) 


n 



w 

bbT 1 


S=1 


< 1 . 


( 6 ) 


Proof. The necessity is simple: for an arbitrary cycle we multiply the n inequalities 
Q-isQisls i < a i s +n s = 1, • ■ ■, n, and obtain (0. To prove sufficiency, we slightly increase all 
numbers q\d so that © still holds for all cycles. This is possible, because the total number of 

cycles is finite. We set ay = 1 and ay = max R [” =1 qfj g x , where maximum is computed aver 
all paths i\ —* • • • —* i n —> i n+ 1 with R = 1, i n+ \ — j, n > 0. Note that if a path contains a 
cycle, then removing it increases the product 11^=1 s i nce the corresponding product 

along the cycle is smaller than one. This means that, in the definition of aq, it suffices to 
take the maximum over all simple (without repeated vertices) paths, i.e., over a finite set. 

It is easy to see that caq^f 1 < oq. Reducing now all back to the original values, we 
obtain strict inequalities. □ 


Theorem 2 If the products IR,... IR are dominant and have unique simple leading eigen¬ 
values, then they have an admissible balancing vector. 

Proof . In view of Lemma [31 it suffices to show that for every cycle (R,..., i n ), n > 2, on the 
set { 1 ,..., r}, we have R [" =1 Qi s is +1 < 1 - We denote this quantity by h and take arbitrary 
5 > 0. There is a product II of matrices from A such that | ( 'v* , ) > q ili2 — 5. Without 

loss of generality we assume that this scalar product is positive. Since the product IR 2 has a 
unique simple leading eigenvalue 1 , it follows that for every igl^we have —$■ (v* 2 ,x) v i2 
as k —> oo. Applying this to the vector x = IL;^, we conclude that Hn^nu^ — q^v^W < 25, 
whenever k is large enough. Thus, for the product II (1) = IT^II, the vector is close 

to q ili2 v i2 . Analogously, for each s — 1,..., n, we find a product IT 3 ) such that the vector 
m s) v is is close to qi 3 i s+1 Vi s+1 . Therefore, for the product B = th e vector Bv tl 

is close to (n:= i Qisis+i) v h = h v n■ Note that ||R|| < M , where M is the supremum of 
norms of all products of matrices from A. If h > 1, then invoking Lemma Q] we conclude 
that p(B) > 1 — e, where £ > 0 can be made arbitrarily small by taking k —>■ oo. Due to 
the dominance assumption, it follows that B is a power of some n 0 G fl, where fl is the set 
of products IIi,..., n r and of its cyclic permutations. Due to the dominance assumption, it 
follows that B is a power of some Ido G D- Taking k large enough we apply Lemma [2] to the 
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words a = II 0 , b = H i2 and conclude that n i2 is a cyclic permutation of n 0 . Similarly, n i3 is 
a cyclic permutation of n 0 . This is impossible, because i 2 ^ 13 , and the products hh 2 , Lh 3 are 
not cyclic permutations of each other. The contradiction proves that h < 1 which completes 
the proof of the theorem. □ 

Remark 3 In a just published paper [25], Moller and Reif present another approach for 
the computation of joint spectral radius. Developing ideas from [23] they come up with an 
elegant branch-and-bound algorithm, which, in contrast to the classical branch-and-bound 
method [16] , can find the exact value. Although its running time is typically bigger than for 
our invariant polytope algorithm HZ] (we compare two algorithms in Example (]2]) below), 
it has several advantages. In particular, it uses the same scheme for the cases of one and of 
many s.m.p. It would be interesting to analyze possible application of the balancing idea for 
that algorithm. 

3.4 How to find the balancing vector 

Thus, an admissible balancing vector a does exist, provided our candidate s.m.p.’s are 
dominant products and their leading eigenvectors are unique and simple. To find a, we 
take some k > 0, compute the values q\S by evaluating polytopes P*^, i = 1 ,...,r, set 

(k) (k) 

Hi = log oti , b\j = — log qlj and solve the following LP problem with variables yo -,..., y r - 
f 2/o max 

l lh - Vj < -Vo + b[f, i,j = 1,..., r, i^j. 

If 2/o < 0, then the fc-admissible vector does not exist. In this case, we have to either 
increase k or find other candidate s.nr.p.’s. If yo > 0, then we have a fc-admissible vector 
a = ( e yi ,..., e Vr ). This vector is optimal in a sense that the minimal ratio between hr anc l 

qff over all i , j is the biggest possible. 

Remark 4 To find an admissible vector one needs to solve LP problem ([7]) for k = 00 . 
However, in this case the evaluation of the coefficients may, a priori, require an infinite 
time. Therefore, we solve this problem for some finite k and then run Algorithm Q] with 
the obtained balancing vector a. If the algorithm terminates within finite time, then a is 
admissible indeed (Theorem [T]). Otherwise, we cannot conclude that there are no admissible 
balancing and that our candidate s.m.p.’s are not dominant. We try to to increase k and 
fold a new vector a. 

Thus, Step 5 of Algorithm |J] consists in choosing a reasonably big k and solving LP 
problem (J7J). If it results yo < 0, then the balancing vector does not exist, and hence the 
algorithm will never converge and we have to find another candidate s.m.p. If yo > 0, then 
the vector a = (e yi ,..., e Vr ) is fc-admissible. If the algorithm does not converge with this a, 
we increase k and solve ([T]) again. 
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Remark 5 Our approach works well also if the family has a unique s.rn.p. IR, but there 
are other simple products IR,..., II r for which the values [p{tlj)} 1 ^ nj although being smaller 
than [p(lR )] 1 /™ 1 = 1, are close to it. In this case the (original) invariant polytope algorithm 
sometimes converges slowly performing many iterations and producing many vertices. This 
is natural, because if, say [p(IR )] 1//n2 = 1 — 5 with very small <R 0 , then the dominance of 
hi! over n 2 plays a role only after many iterations. Our approach suggests to collect all 
those “almost s.rn.p. candidates” IR,..., IR add them to IR, find the balancing multipliers 
i a i}i=i f° r their roots by solving LP problem ([7} and run Algorithm [T] for the initial set 
Vo = {aj'Hj} r j =1 . In most of practical cases, this modification significantly speeds up the 
algorithm. 

Another modification of the invariant polytope algorithm is considered in the next section. 

Example 2 Consider the following example introduced by Deslaurier and Dubuc in (E|, 
associated to an eight-point subdivision scheme, 
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The joint spectral radius of A = {A 4 ,A 2 } was found in [22], where it was shown that 
both A\ and A 2 are s.rn.p. Its computation required the construction of a binary tree with 
14 levels and considering about 130 matrix products (i.e., vertices of the tree). Applying our 
Algorithm Q] with the candidates s.rn.p. Ai and A 2 and with a balancing vector a = (1 1 ) T 

for the leading eigenvectors V\ and v 2 of A 4 and A 2 , respectively, we construct the invariant 
polytope with 24 vertices in 5 steps: 


V\ 


V2 


V 3 = 

Aiv 2 

v 4 = 

= A 2 V i 

V 5 -- 

= Ai v 3 

V6 = 

= A, v 4 

v 7 -- 

- A 2 v 8 

V 8 ~- 

- A 2 v 4 

v 9 = 

-- Ai v 5 

Vio 

= A.i v 6 

vn 

— Ai v 7 

V 12 

= A\ v 8 

V 13 

— A 2 u 5 

V\4 

— A 2 Vq 

V 15 z 

= A 2 v 7 

Vie 

— A 2 v 8 

Vyj 

— Ai vn 

Vl8 

— A-i v 12 

Vl9 

= Ai Vi3 

V20 

= Ai Vi4 

V21 : 

= A 2 v n 

V22 

= A 2 V\2 

V 23 

= A-2 U 13 

V24 

= A 2 v u 
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Thus, in our case it suffices to construct a binary tree with 5 levels and consider 24 of its 
vertices. 


4 Introducing extra initial vertices 

The same approach developed for the case of many s.m.p. can be used to introduce extra 
initial vertices. Sometimes Algorithm Q] converges slowly because the family A is not well- 
conditioned: its matrices have a common “almost invariant subspace” of some dimension s < 
d — 1. In this case the invariant polytope P may be very flattened (almost contained in 
that subspace). As a consequence, the algorithm performs many iterations because the 
polytopes Pfc, being all flattened, badly absorb new vertices. To avoid this trouble one can 
introduce extra initial vertices xi,...,x s and run Algorithm Q] with the initial set Vo = 
{a{Hi ,..., a r H r , xi ,..., x s }- In next theorem we use the value qj(x) defined in (HI) . 

Theorem 3 Suppose Algorithm^ with initial roots Hi ,..., Hr terminates within finite time; 
then this algorithm with extra initial vertices x ±,..., x s also does if and only if qfixf) < 1 for 
all j — 1 ,..., r; i — 1 ,..., s. 

Proof . For the sake of simplicity, we consider the case s = 1, r = 1, the proof in the general 
case is similar. Let P denote the final polytope produced by the algorithm starting with the 
root Hi, and let P(xi) = absco {flxi | fl e A n , n > 0}. 

Necessity. Assume the algorithm terminates within finite time. In the proof of Theorem Q] 
we showed that the maximum of the linear functional f(x) = (u*,x) on the final polytope P 
is equal to one and is attained at a unique point Vi. Hence, either qfixi) < 1, in which case 
the proof is completed, or qfixfi = 1, and hence there is a sequence of products {IT (fc ^} fceN 
such that K,IT%i) —> 1 as k —* oo. Therefore, IV fe 'xi —* Vi as k —* oo. This implies that, 
for sufficiently large k, the points IV^Xi are not absorbed in the algorithm, and hence, the 
algorithm cannot terminate within finite time. 

Sufficiency. Since the second largest eigenvalue of the matrix Hi is smaller than 1 in 
absolute value, it follows that n^' —$■ V\ [u^]- 71 as k —> oo. If qfixfi < 1, then the matrix 
Vi [vf] 1 maps the set P(x i) to the segment [— qi(xi)vi, qi(xi)vi], which is contained in 
qfixfiP. Hence, n^' (P(xi)) C P, for some k. Therefore, every product n of length N 
containing a subword takes the point xi inside P. On the other hand, for all products 
n not containing this subword, we have ||n|| < Cq N , where C > 0, q G (0,1) are some 
constants (see m Theorem 4]). Hence, for large N, all such products also take the point Xi 
inside P. Thus, all long products take Xi inside P , hence the algorithm starting with the 
initial set Hi U {aq} terminates within finite time. □ 

In practice, it suffices to introduce extra initial vertices of the form x* = /^e*, where e; 
is the canonical basis vector (the i-th entry is one and all others are zeros) and /q is some 
positive coefficient. We fix a reasonably small £ > 0, say, between 0.001 and 0.1, a reasonably 
large k, say k = 15, rum k iterations of Algorithm |Tj and compute the values 

Qi = max I (e*, v) I , i = 1 ,..., d . 
vev k 1 v 5/1 
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(k) 

So, Q\ is the length of projection of the polytope onto the i-th coordinate axis, or the 
largest i-th coordinate of its vertices. If Q\ > e for all i, then P% contains the cross-polytope 
absco { eei , i — 1,..., 5}, which, in turn, contains the Euclidean ball of radius ej\fd centered 
at the origin. In this case the polytope P *. is considered to be well-conditioned, and hence 

(k) 

we do not add any extra vertex. If, otherwise, Q\ < e for some i, then we add an extra 
vertex Xi = e [max.,- g^(ej)] -1 ej. Collecting all such vertices for % — 1,..., d (assume there 
are s < d ones) we run Algorithm Q] with s initial extra vertices in the set Vo- 

In Section Owe apply this trick to speed up Algorithm [Q for Daubechies matrices, which 
turn out to be extremely ill-conditioned. 

5 Applications: the Butterfly subdivision scheme 

Subdivision schemes are iterative algorithms of linear interpolation and approximation of 
multivariate functions and of generating curves and surfaces. Due to their remarkable prop¬ 
erties they are widely implemented and studied in an extensive literature. 

The Butterfly scheme originated with Dyn, Gregory, and Levin [2] and became one of 
the most popular bivariate schemes for interpolation and for generating smooth surfaces 
(see also the generalization given in [38]). This scheme is a generalization of the univariate 
four-point interpolatory scheme to bivariate functions [HE31E5]. First we take an arbitrary 
triangulation of the approximated surface and consider the corresponding piecewise-linear 
interpolation. This interpolation produces a sequence of piecewise-linear surfaces with thri- 
angular faces that converges to a continuous surface, which is considered as an interpolation 
of the original one. To describe this algorithm in more detail we assume that the original 
surface is given by a bivariate function f(x 1,0:2). We consider a regular triangulation of R 2 
and take the values of the function / at its vertices. So we obtain a function j\ defined on 
a triangular mesh. In the next iteration we define the function / 2 on the refined triangular 
mesh: the values at the vertices of the original mesh stay the same, the values at midpoints 
of edges are defined as a linear combination of eighth neighboring vertices as shown in the 
following figure (X is the new vertex, the coefficients of the linear combination are written 
in the corresponding vertices). The parameter oj is the same for all vertices and for all itera¬ 
tions. In the next iteration we do the same with the new mesh and produce the function / 3 , 
etc. Each function is extended from the corresponding mesh to the whole plane by linear¬ 
ity. The scheme is said to converge if for every initial function /, the functions fk converge 
uniformly to a continuous function 5(/). Due to linearity and shift-invariance, it suffices 
to have the convergence for the initial 5-function, which vanishes at all vertices but one, 
where it is equal to one. The corresponding limit function 5(5) is called reftnable function 
and denoted by Lp. The scheme converges for every uj 6 (0, ^) and reproduces polynomial 
of degree one, i.e., if f(x i,x 2 ) = cqaq + a 2 x 2 + ao, then 5(/) = /. The case c 0 = 1/16 is 
special, in this case the scheme reproduces polynomials of degree 3 na. One of the most 
important problems in the study of any subdivision scheme is its regularity. We use the 
standard modulus of continuity 

ujf(h) = sup{\f(x + f) - f(x)\ | *eRM|£||<h}. 
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Figure 1: The generalized Butterfly scheme. 

The Holder exponent of the function / is 

ctf(h ) = n + sup{a > 0 | cofoo (h) < Ch a , h > 0 } , 

where n is the biggest integer such that / E C n (R d ). The Holder regularity of a subdivision 
scheme is a v . This is well-known that the exponent of Holder regularity of a bivariate 
subdivision scheme is equal to 

a* = -log 2 p(Tf,T 2 W ,Tf,Tf), 

where £ is the maximal degree of the space of algebraic polynomials reproduced by the 
scheme, if' 1 are restrictions of the transition operators Ti of the scheme to their common 
invariant subspace orthogonal to the subspace of algebraic polynomials of degree i [.26]. For 
the Butterfly scheme, for all u> 1/16, we have £ = 1 and the operators Tf' are given 
by 24 x 24-matrices. In the only “most regular” case oj = 1/16 the scheme respects the 
cubic polynomials, i.e., i = 3, and T> ’ are 17 x 17-matrices. For this exceptional and most 
important case it was conjectured in early 90th that a^ = 2, i.e., the limit function p of the 
scheme are continuously differentiable and its derivative p 1 has Holder exponent 1. We are 
going to prove this conjecture and, moreover, we show that the derivative of the function p 
is not Lipschitz, but “almost Lipschitz” with the logarithmic factor 2. 

Theorem 4 The Holder regularity of the Butterfly scheme with cu = -f is equal to 2. The 
derivative p l of the limit function is “almost Lipschitz’’ with the logarithmic factor 2: 

c<V(h) x h | log h\ 2 , he(o,^). (8) 

Remark 6 In the four-point subdivision scheme, which is a univariate parameter-dependent 
analogue of the butterfly scheme, the case uj = ^ is also crucial. This is the only case when 
the scheme reproduces cubic polynomials. As it was proved by S.Dubuc in 1986 jl3j . the 
regularity of the four-point scheme in this case is equal to two and co^flh) x h\ logh|, i.e., p' 
is almost Lipschitz with the logarithmic factor 1. By Theorem [4] for the Butterfly scheme, 
the situation is similar, but p' is almost Lipschitz with the logarithmic factor 2. 
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To prove Theorem [4] we first show that p{T^ 3 \ i = 1,..., 4} = |. Then we conclude that 
p{T) , i — 1,..., 4} = /. By a more refine analysis of the matrices we establish that 


max 


7^(1). # # 7 ^(i) || 
± i k ± h I 


h, ■ ■ ■ ,h £ { 1 , 2 , 3 , 4 } | 


k 2 4~ k . 


( 9 ) 


Then it will remain to refer to some known facts of the theory of subdivision schemes. The 
main and most difficult part is the finding of the joint spectral radius of the matrices T- 3 \ 
This is done in the next subsections. Then we conclude the proof. 


5.1 The case uj = 1/16: the classical Butterfly scheme 

The 17 x 17-matrices T® can be computed exactly by the Matlab program of P.Oswald [27]. 

/o\ 

To simplify the notation, we denote Ai = 4 T> , i = 1 ,..., 4, and A = {7b, ..., A 4 }. These 
four 17 x 17-matrices are written in Appendix 2. Our goal is to prove that p(A) = 1. It 

is remarkable that all the matrices of the family A and the leading eigenvector of its s.rn.p. 

possess rational entries, so our computations are actually done in the exact arithmetics. 

Step 1. Factorization of the family A to A\ and „4. 2 . 

It appears that the matrices Ai ,..., A 4 can be factored to a block lower-triangular form 
in a common basis. To see this we take the leading eigenvectors V\ of the matrix A\ (as it 
was mentioned above, it has rational entries): 

_ ( 1746890 942260 1391945 1596995 

11 2004757 ’ 2004757 ’ 2004757 ’ 2004757 ’ 

1478789 2902096 594645 1063787 2802569 2242099 
2004757 ’ 2004757 ’ 4009514 ’ 4009514 ’ 4009514 ’ 4009514 ’ 

It is checked directly that the linear span of the following six vectors is a common invariant 
subspace of all matrices from A: Vi, v 2 = A 2 v i, v 3 = A 3 v i, n 4 = Api!, v 5 = A 1 A 2 v 1 , v 6 = 
AiA 3 v\. Therefore, we can transform the matrices from A into block lower-triangular form 
with diagonal blocks of dimensions 6 and 11. The transformation matrix is 

S = { e i e 2 ... en v x v 2 v 3 v 4 v 5 v 6 ) 

(written by columns), where e*, G M 1 ' is the k -th canonical basis vector. This gives the 
transformed matrices with block lower-triangular structure: 

S 1 A i S =( D . B . ) ’ ? = 1,...,4. 

with 6x 6-matrices C) and 11 x 11-matrices Bi. Those matrices are written down in Appendix 
2. Note that they are all rational. Let ^4i = {£>i, B 2} B 3 , £> 4 } and ^4 2 = {Ci, C 2 , C 3 , C* 4 }. It 
is well known that the joint spectral radius of a block lower-triangular family of matrices is 
equal to the maximal joint spectral radius of blocks [3]. Hence p(A) = max {p(Ai), p(^4 2 )}. 
Now we are going to show that p(A\) = p(A 2 ) = 1, from which it will follow that p(A) = 1. 
We begin with the family A \. 


1987045 

’4009514 ’ 

45664 

'2004757 ’ 


3186205 1392363 


4009514 ’ 4009514 ’ 

o.o, A 
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Step 2. Analysis of the family A\ 

We have a family A\ of 11 x 11-matrices , B 4 written in Appendix 2. Each of the 

matrices B 1, B 2l B 3 has a simple leading eigenvalue 1, the corresponding leading eigenvectors 
Mi,m 2 ,m 3 are all simple. The matrix B A has spectral radius 1/2. We are going to show that 
p(A\) = 1, i.e., this family has three s.m.p.: B\, B 2 and B 3 . The leading eigenvectors are 
(normalized in the maximum norm), 


Ml 


2497 \ 
3306 
1453 
3306 
1997 
3306 
283 
3306 
2023 
3306 
787 
3306 
2851 

11020 

4749 

11020 


U 2 


( 


230 \ 

1089 ' 

230 
1089 
370 

1089 

730 

1089 

730 

1089 

370 

1089 


U3 


( 


66 \ 
625 
1158 
4375 
2298 
4375 
82 

4375 

154 

625 

442 

4375 


221 

363 

1 


-1 

2661 

21875 


613 

2204 

-1 


V 1 


221 

363 

0 


15959 

21875 

2204 

4375 


V 0/ 


2204 
4375 / 


Solving problem (J7J) for r = 3, k 


10, and for the family Ai, we obtain the scaling factors: 


a\ = 0.50379 ..., a 2 = 0.48126 ..., 0:3 = 1 


The Algorithm Q] with the three candidate s.m.p.’s IIi = Bi,U 2 = B 2 , and II3 = B 3 and 
with those factors a* terminates within four iterations and produces an invariant polytope. 
However, we slightly change the factors to a\ = 0.5, a 2 = 0.5 and a .3 = 1 respectively 
in order to preserve the rationality of the vectors (and consequently the exactness of the 
computation). Thus, 

1 1 

V\ = “Mi, V 2 = —M 2 , M 3 = M 3 . 

The algorithm still converges within four iterations producing the polytope P\ with 75 • 2 
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vertices. Here is the list of vertices: 


V\ 


V2 


V 3 


V4 -- 

= Bi v 2 

V 5 -- 

= Bx v 3 

V6 z 

= B 2 Vi 

V 7 -- 

= B 2 v 3 

Vs -- 

= B 3 v 1 

Vg ~- 

= B 3 V 2 

VlQ 

= b 4 

Vx 

Vll 

= B 4 v 2 

V12 

= B 4 v 3 

VX 3 

= B1V4 

Vl 4 

= Bx v 5 

V15 

= Bx 

v e 

vie 

= B 1 v 7 

vn 

= Bx v 8 

Vis 

= Bx Vg 

Vig 

= Bx v 10 

V20 

= Bx 

vn 

V21 

= V12 

V22 

= B 2 V 4 

V23 

= B 2 v 5 

V 2 A 

= B 2 v 6 

V 25 

= b 2 

v 7 

V26 

= B 2 v 8 

V27 

= B 2 Vg 

V 28 

= B 2 v 10 

V 2 9 

= B 2 vn 

V30 

= b 2 

V12 

V31 

= B 3 v 4 

V32 

= B 3 v 5 

v 33 

= B 3 v e 

V34 

= B 3 V7 

V 35 

= b 3 

Vs 

v 3 e 

= B 3 Vg 

V37 

= B 3 vig 

V 3 8 

= B 3 vn 

v 3 g 

= B 3 VX2 

V 40 

= b 4 

V4 

V41 

= B 4U5 

V42 

= B 4 v 6 

V43 

— B4 V7 

V44 

= B 4 v 8 

V45 

= b 4 

Vg 

v 4 6 

— B4V1Q 

V47 

= B4 v u 

V48 

— B 4 Vi 2 

V49 

— Bx Vx 5 

v 50 

= Bx 

Vl 7 

V51 

= Bi v 20 

V52 

= Bx v 2 1 

V 53 

= Bx v 30 

V54 

= Bx v 38 

V55 

= Bx 

V46 

V56 

= B 2 v 2 \ 

V57 

= B 2 v 22 

V 5 8 

= B 2 V 2 7 

V59 

= B 2 v 28 

veo 

= b 2 

v 3 o 

vei 

= B 2 V 3 j 

VQ 2 

= B 2 V47 

V &3 

= B 3 V 20 

VQ 4 

= B 3 v 28 

V &5 

= b 3 

V 3 2 

vm 

— B 3 V 3 4 

VQ 7 

— B 3 v 37 

V 6 8 

= B 3 v 38 

Vgg 

= B 3 V 48 

v 70 

= B 4 

V 20 

V71 

— B 4 v 2 \ 

V 7 2 

= B 4 v 28 

v 73 

= B 4 V 30 

V 7 4 

— B4 v 37 

V 75 

= b 4 

V 3 8 


Thus, p(Ai) = 1 and, by Theorem [lj B ll B 2l B 3 are dominant products for Ax- 

Step 3 . Analysis of the family A 2 

We have a family A 2 of 6 x 6 -matrices C'i,..., C 4 written in Appendix 2 . Each of the 
matrices C±, C 2 , C 3 has a simple leading eigenvalue 1 , the matrix C4 has spectral radius 1 / 2 . 

First of all, we observe the existence of three invariant 2 -dimensional subspaces of all the 
matrices C', : . We indicate by w 7 , w 2 and w 3 the unique leading eigenvectors associated to the 
eigenvalue 1 of the matrices C\C 2 , C\C 3 and C 2 C 3 , (normalized in maximum norm), 


wi = 


/ i \ 

1 

4 

0 
0 

V 1 / 


w 2 = 


(\\ 


1 

4 

0 

0 

V 1 / 


-n 


w 3 = 


25 

'28 

25 

'28 

_ 2 
7 


V 


The invariant subspaces are given by V\ = span (uq, C4W1), V 2 = span (w 2 , C 4 w 2 ) and V 3 = 
span(w 3 ,C 1 w 3 ). Thus we define the matrix 


S = (wi, C4W1, w 2 , C 4 w 2 , w 3 , C\w 3 ) 


which block-diagonalizes all matrices Cj, i = 1 ,... , 4 . We denote the diagonal blocks of the 
matrices S~ 1 CjS as G ,,\, G ; 2 and GA. We obtain three families of 2 x 2 matrices to analyze, 
G^iG^G^G^Gu} with 
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then G2 — {(?2i, G22, G23, G24} with 



i.e. G2 — Gi and G3 — {G31, G32, G33, G34} with 



All previous families have joint spectral radius 1. The Li-norm is extremal for G \. This 
means that ||Gij||i < 1 for j — 1 ,..., 4. Hence p(G\) = 1. Since G 2 = Gi, it follows that 
p((j2) = 1- For the family G3, we apply Algorithm [I] and obtain the invariant polytope P. 
In this case P is an octagon with vertices 



Thus, p(Gz) = 1, and hence p(A.2) = max{p(Gi),p(G2),p(G3)} = 1. 

The proof of Theorem |4| 

We start with introducing some further notation. A norm || • || in W l is called extremal 
for a family A if ||Aj|| < p(A) for all Aj e A. Algorithm [Q constructs an extremal polytope 
norm. A family A is called product bounded if norms of all products of matrices from A are 
uniformly bounded (see e.g. [2T]). If a family has an extremal norm and p(A) = 1, then it 
is product bounded. 

We have shown that p(A\) = p(A 2 ) = 1. Hence, the block lower-triangular form yields 
that p(A) = max{p(Ai), p(A.2)} = 1, and so p{T®, i — 1,..., 4} = \ . Furthermore, all 
matrices Tw in a special basis of the space M 21 have the form: 

/ J 2 0 0 \ 

T. (1) = * J 3 0 , <=1,2,3,4, (10) 

V * * G ) 

where J s is the (s + 1) x (s + l)-diagonal matrix with all diagonal entries equal to 2~ s 
(see HD 1ZS3)- Therefore, the joint spectral radius of {T^G = 1, ...,4} is equal to the 
maximum of the joint spectral radii of the three blocks, i.e., the maximum of p(J 2 ) = 

of p(Js) = |, and of p{T,- 3 \i = 1,... ,4} = j. Thus, p{T^\i = 1,... ,4} = and hence 
ay = — log 2 | = 2. The Holder exponent is found. Now let us analyze the regularity of the 
derivative <p'. 

For any refinable function <p, the modulus of continuity cu^i^h) is asymptotically equiv¬ 
alent to the logarithm of the left-hand side of the equality (J2J) with k = — [log 2 h] (see |31j). 
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Hence, to prove that (J v >(h) x h\logh\ 2 it suffices to establish (ED- Applying factoriza¬ 
tion (HU and the results of Steps 1-3, we obtain 



/ J 2 0 0 0 \ 

* ,/y 0 0 

* * \Bi 0 

V * * * i c i ) 


i = 1,2, 3,4. 


In this block lower-triangular form, we have three blocks (J 2 , \Bi and \Cj) with the joint 
spectral radius 4 and one (J 3 ) with a smaller spectral radius (|). Moreover, all these former 
three blocks are product bounded, since they have extremal norms. Therefore [31] . 

max |||| < C^k 2 , ke N , (11) 


where C\ is a constant. On the other hand, it is verified directly that each of the matrices 
T- l \ i = 1, 2, 3, has two Jordan blocks of size 3 corresponding to the leading eigenvalue A = 
Hence, the left-hand side of (HIT) is bigger than or equal to || [Xj (1) ] fc || > C 2 X k k 2 = 4 ~ k k 2 . 
Therefore, it is asymptotically equivalent to 4 ~ k k 2 . This proves (ED and hence x 

h|logh| 2 . □ 


5.2 Other values of the parameter uj 

The convergence analysis of the Butterfly scheme can be extended to other values of uj G 
[0,4]. In this case we have to deal with 24 x 24-matrices T^\i = 1,2, 3,4. The scheme 
converges (to continuous limit functions) if and only if their joint spectral radius p is smaller 
than one. The regularity of the scheme is equal to = — log 2 p. 

For 0; = ^ each matrix T^\i = 1,2,3, has two simple eigenvalues of modulus one: 
precisely 1 and —1. The matrix T 4 ' 11 has a simple eigenvalue —1 and the 1 of multiplicity 2 
(both algebraic and geometric). The two leading eigenvectors corresponding to 1 and —1 
define a common invariant subspace for the family which can be transformed into a similar 
block triangular form. The 2 x 2-blocks are respectively 



They are all symmetric, hence their joint spectral radius equals to the maximal spectral 
radius of these matrices [9], i.e., is equal to one. The remaining 22 x 22 family of 4 matrices 
has the fourth matrix as an s.m.p. Starting from its (unique) leading eigenvector Algorithm|T] 
terminates within 8 iterations and constructs an invariant polytope norm with 487 vertices. 
This proves that p(T^ ,..., T 4 1 ) = 1, and hence the scheme does not converge. 

We have successfully applied our procedure also for other values to G [0, |). This leads 
us to conjecture that the generalized Butterfly subdivision scheme is convergent in the whole 
interval. To support this conjecture we report in Table Q] the results obtained for c 0 = 
k = 0,1,..., 15,16 (see also Figure ED- 
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k 

0 

1 

2 

3 

4 

5 

6 

7 

8 


1 

1.1000 

1.2284 

1.4150 

2 

1.6781 

1.4150 

1.1926 

1 

k 

9 

10 

11 

12 

13 

14 

15 

16 



0.8301 

0.6781 

0.5406 

0.4150 

0.2996 

0.1926 

0.0931 

O' 



Table 1: Computed Holder exponent of Butterfly scheme for uj k = k = 0,..., 16. 



Figure 2: The computed Holder exponent of the generalized Butterfly scheme. 

6 Applications: the regularity of Daubechies wavelets 

One of the most important applications of the joint spectral radius is the computation of the 
Holder regularity of refinable functions and wavelets. For Daubechies wavelets, this problem 
was studied in many works (see 0 ei mu hd ng E0 S3 ezi and references therein). Let 
us recall that the Daubechies wavelets is a system of functions 2^ 2 ip(f2^x — n), j,n G Z, 
that constitutes an orthonormal basis in All functions of this system are generated 

by double dilates and integer translates of the compactly supported wavelet function if. 
I.Daubechies in pm constructed a countable family of wavelet functions ip — 'f’N, N > 1, 
each generates its own wavelet system. The function ipi is the Haar function. For all N > 2 
the functions ipN are continuous, their smoothness increases in N and a^ N > 0.2 N [TO]. So, 
there are arbitrarily smooth systems of wavelets. However, the price for the regularity is the 
length of the support, which also grows with N\ supply = [0, 2 N — 1], The regularity is a 
very important characteristics of wavelets, in particular, for their applications in functional 
analysis, approximation theory, image processing and in numerical PDE. There are several 
methods to obtain lower and upper bounds for the Holder exponents of wavelet functions 
(see 0 HDl ESI E31E3). The matrix approach is the only one that theoretically allows to 
find them precisely. It was established in U HU that a^ N = N — log 2 p(H 0 , Hi), where 
Ho, Hi are special matrices of size (N — 1) x (N — 1). This enabled to find the precise 
values of the Holder exponent for some small N. For N = 2, 3, and 4, the value a^ N were 
found by Daubechies and Lagarias in HP; for N = 5,6, 7, and 8, they were computed by 
G.Gripenberg [T6]. Every time a delicate analysis of special properties of those matrices was 
involved. In all the cases the s.m.p. of the family {Ho, Hi} was one of those two matrices, 
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and it was a general belief that this is the case for all N. We apply the standard routine of 
Algorithm Q] to find the precise values of a^ N for all N < 20. In particular, we shall see that 
for N = 10, the conjecture of one matrix s.rn.p. is violated and the s.rn.p. is BqB\. 

We need to recall key steps of construction of the matrices Bq,B\. For every N = 1, 2,... 
we have a set of 2N Daubechies filter coefficients : Co,..., C2jv-i- They possess some special 
properties, in particular, ^^ 0 _1 c* = 2 and the polynomial m(z) = 1 CnZU ^ as zero °f 

order N at the point z — — 1. We set 


q(z) 


m(z) 

((1 + Z )/ 2 ) N 


N -1 


n=0 


and write the transition k x /e-matrices as follows: 


{Boh = Qx-j-u ( B ih = 92i-j, i, j = 1, ■ ■ •, N - 1. 

We compute p(Bo, Bf) by Algorithm [Q For some N we have a non-unique s.rn.p. (these are 
the cases when Bq and B\ are both s.rn.p.) and find the balancing vector a by the method 
in Subsection 13.41 However, due to symmetry, the entries of that two-dimensional vector a 
are equal. Another difficulty, much more significant for the Daubechies matrices is that 
with growing N they become very ill-conditioned. All vertices of the constructed polytope 
have very small last components, which corresponds to the property of quasi-invariance of a 
certain subspace and determines a polytope strongly flattened along certain directions. For 
N = 10, the last components are about 10 -12 — 10~ 19 of the values of the first components. 
This creates enormous numerical difficulties in the running of Algorithm [T] in particular, in 
the linear programming routines. That is why we use the technique with extra initial vertices 
(Section H]). In the next subsections we present three illustrative cases (N = 4,10,12) and 
report the computed Hdlder regularity of Daubechies wavelets for al IV < 20. 


6.1 Illustrative examples 

We demonstrate the computation process for IV = 4,10 and 12 and see the crucial changes 
in the behaviour of Algorithm [T| when the dimension grows. The case N = 4 was done (by a 
different approach) by Daubechies and Lagarias in [11], while the two other cases are new. 

The case N = 4. 

For the pair of transition matrices: 

/ 5.212854848820774 0 

B 0 = 1.703224934278843 -4.676287953813834 

\ 0 -0.239791829285782 

/ -4.676287953813834 5.212854848820774 

B x = -0.239791829285782 1.703224934278843 

\ 0 0 


°\ 

5.212854848820774 , 

1.703224934278843 / 

°\ 

-4.676287953813834 , 

-0.239791829285782 / 
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the candidate s.m.p. is B 0 with p(B 0 ) = 5.212854848820774 .... In order to apply Algorithm 
[Owe compute the leading eigenvector of B 0 , 


/ 1.0000 \ 

0.1662 

\ -0.0113 / 


and set Vo = {r>i}. Observe that v\ almost lies on the subspace E 2 C R 3 spanned by the 
vectors {ei,e 2} of the canonical basis of M 3 . Applying the normalized matrices B 0 and B\ 
repeatedly to Vo one observes that the resulting vectors also almost lie on the subspace E 2 . 
This has implications on the flatness of the invariant polyhedron computed by Algorithm 
Q] In its basic implementation the algorithm converges and generates a centrally symmetric 



Figure 3: The invariant polytope for the Daubechies matrices for N = 4 computed by the 
standard algorithm. 


polytope in M 3 of 6 • 2 vertices, in 7 iterations. The partial polytope norms of the family 
A = {B 0 , Bi} are reported in Tabled The vertices (beyond V\) of the polytope follow (we 


Iteration k 

4 5 6 7 

1 ' II Vi 

1.7845 1.1288 1.0571 1 


Table 2: The partial polytope norms computed by the standard algorithm for k — 4. 


0.7308 \ 
0.0108 , 
0.0115 / 

/ -0.7308 \ / -0.7308 \ 

-0.2175 , v 6 = -0.0393 

\ 0.0042 / \ 0.0113 / 
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report only half of them): 

-0.7308 
v 2 = | 0.0185 


5.2249 • 10 


-4 


V 3 = 


0.7308 
0.2548 
6.8062 • 10" 4 


v 4 = 


v 5 = 












and the corresponding unit polytope P = abscojwj, i — 1,..., 6} is shown in Figure [3] We 
see that P appears to be very flat. 



Figure 4: Polytope extremal norm for the Daubechies matrices for N = 4 computed with an 
extra initial vertex. 

In fact, the largest singular valuc@ of the matrix of its vertices V = {±fj, i = 1,..., 6} is 
cri = 1.9365, while the smallest one is <73 = 1.0923 • 10” 1 2 (note that if a 3 would be zero then 
V would not span the whole space and the polytope would be contained in a subspace). This 
almost 200 times difference gives a numerical evidence of the flattening phenomenon. To 
avoid it we add an extra initial vertex x\ = ^e 3 and obtain a better behaviour of Algorithm 
[Hand a more balanced polytope (see Figure [4]). The polytope has now 4 • 2 vertices which 
are computed in only two iterations. Denote v 2 = X\. The vertices beyond v\ and V2 of the 
polytope follow (we report only half of them): 

( °\ 

v' 3 = 0.8000 , v' A = 

\ 0.2613 / 

The largest singular value is now o\ = 1.1597 • 10° and the smallest singular value is now 
<j 3 = 7.7402 • 10 -1 , which demonstrate a much more balanced shape of the unit ball of the 
polytope extremal norm. The computed Holder exponent of is = 4 — log 2 p(B 0 , B\) = 
4 - log 2 p(B 0 ) = 1.6179.... 

The case N = 10. 

We have the 9 x 9-matrices B 0 , B 3 and are going to prove that the s.m.p. is 
hi = Bl Bf, p c = p(n ) 1/4 = 99.636965469277555 ... 

1 Recall that for a matrix B £ (or B £ C p,q ), the reduced singular value decomposition is given by 

B = UT,W* where U £ C p,q and W £ C q ’ q are unitary matrices and £ £ R 9,9 is a diagonal matrix with 
nonnegative diagonal elements {tJi}® =1 (the singular values, usually ordered in decreasing way) that are the 
square roots of the eigenvalues of the Hermitian (semi)-positive definite matrix B*B (see e.g. [2011. 


( 0 

-0.7176 

\ -0.0368 
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which interesting in itself since it contradicts to the conjectured property that p(B 0 ,Bi) = 
ma x{p(B 0 ), p(Bi)} for all N (see Introduction). 

Let B 0 = B 0 /p c , Bi = Bi/p c . Applying Algorithm [Q we compute the starting set of 
vectors Vo = {wi, v 2 , V3, rq}, where v\ is the leading eigenvector of II, v 2 = Bpvi, V3 = 
BiV 2 ,V 4 = B 0 V3. Clearly, B0V4 = v\. Thus, v 2 ,V3 and V4 are the leading eigenvectors of 
cyclic permutations of the product II. We have 


1.7122 • 10" 3 \ 


/ 3.8568 • 10” 1 \ 

1.0000 • 10° 


1.1649 • 10° 

3.0340 • 10” 1 


3.6053 • 10" 2 

-3.0515 • 10- 1 


-1.8840 ■ 10' 1 

-5.9219 • 10" 2 

, v 2 = 

-1.8276 ■ 10" 2 

5.9518 • 10~ 4 


5.0899 • 10" 4 

9.0971 • 10~ 5 


1.0986 • 10" 5 

-2.5228 • 10" 7 


-1.8367- 10" 8 

-6.6280 • 10" 12 ) 


\ 1.2778 • 10~ 15 / 


/ 1.1395 • 10~ 2 \ 


/ 4.4172 • 10" 3 \ 

1.2195 • 10° 


-1.1524 • 10° 

5.6969 • HT 1 


-8.3168- lO' 1 

—2.4167 • 10 -2 


4.6863 • 10" 2 

— 1.2787 • 10 -2 

, V4 = 

3.7374 ■ 10" 2 

-9.4332 • 10” 5 


7.1841 • 10" 4 

4.9101 • 10~ 6 


-3.4331 • 10" 5 

-2.1784 • 10" 9 


3.4387- 10~ 8 

\ -2.4634 • 10" 19 j 


\ 4.1997- 10~ 13 / 


Observe that the last two components of all the vectors v±,... ,V4 are very small, i.e., all 
these vectors almost lie on the subspace E 7 C R 9 spanned by the vectors {ei,..., e 7 } of the 
canonical basis of R 9 . Applying B 0 and Bi repeatedly to Vo one observes that the resulting 
vectors also almost lie on the subspace E 7 . This means that the invariant polytope computed 
by Algorithm Q] is nearly degenerate, it is close to a 7-dimensional polytope. A consequence 
of this is a slow convergence behaviour of the algorithm and an ill-conditioning of basic 
linear algebra operations. The algorithm terminates and generates a centrally symmetric 
polytope in R 9 of 220 • 2 vertices, in 16 iterations. The partial polytope norms of the family 
A = {/!()• B \} are reported in Table El The largest singular value of the set of vertices is 


k 

5 

6 

7 

8 

9 

10 

11 

12 

13 

14 

15 

16 

II ' Ita 

24.856 

4.693 

2.990 

2.237 

1.743 

1.414 

1.140 

1.064 

1.027 

1.025 

1.001 

1 


Table 3: The partial polytope norms computed by the Algorithm Q] for N = 10. 


(T 1 = 1.3253-10 1 and the smallest singular value is <79 = 1.0620-10 10 , which gives a numerical 
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evidence of the flatness of the polytope. To correct the behavior of the algorithm we add 
an extra vector along the “most narrow” (for the polytope P ) direction (see Section [I]). 
Adding the vector n 5 = |eg, we obtain a better behaviour of the algorithm and a more 
balanced polytope. The results are summarized here; the polytope has 75 • 2 vertices which 
are computed in 10 iterations. The largest singular value is now oq = 7.7943 • 10° and the 


i 

2 

3 

4 

5 

6 

7 

8 

9 

10 

1 ' 1 \Vi 

9.4019 

1.8655 

1.3861 

1.2622 

1.2076 

1.1150 

1.0512 

1.0169 

1 


Table 4: The partial polytope norms computed by the modified algorithm for N = 10. 


smallest one is <7g = 1.1401 • 10 3 , which demonstrate a much more balanced shape of the 
unit ball of the polytope extremal norm. The computed Holder exponent of ip io is 

a vi0 = 10 - log 2 p(B 0 , Bi) = 10 - ^ log 2 p{BlBl) = 3.361390821401114... 


Remark 7 If we consider the adjoint family {Bq, B\} we have naturally p(B 0 , Bi) = p(Bq, B^). 
It is remarkable that applying Algorithm Q] to {Bq,BI}, we do not have problems with flat¬ 
tening. Alas convergence remains slow (12 iterations and 370 • 2 vertices). 


The case N = 12. 


In this case both B 0 and Bi are s.m.p., i.e. p(B 0 ,Bi) = p(B 0 ) = p(B i). The balancing 
technique (Section [3]) gives a = (1,1)) he., equal weights to the leading eigenvectors of B 0 
and of Bi. The two vectors, say V\ and v 2 follow: 


0 ^ 


/ -4.8598-10" 1 

1.3465 • 10" 1 


-3.2838 • IO' 2 

1.0000 • 10° 


1.0000 • 10° 

4.3937- 10" 1 


3.5198- 10' 1 

—1.1888 • 10 _1 


-8.3828 • 10~ 3 

-4.1549 • 10" 2 

, v 2 = 

-5.6206 • 10 -3 

-5.5942 • 10" 4 


-5.3176 • IO" 5 

1.2299 • 10" 4 


3.0075 • IO" 6 

2.0148 • IO' 7 


-3.7539 • IO" 9 

-3.6435 • 10~ 9 


-1.2233 • 10~ 13 

1.1240 • IO" 13 j 


V o 


We observe that they almost lie on the subspace E 8 C M 11 spanned by the vectors {e \,..., e 8 } 
of the canonical basis of M 11 . Applying B 0 and Bi repeatedly to Vo one observes that the 
resulting vectors also almost lie on the subspace E 8 . Note that the last components seems 
to vanish exponentially in the number of iterations (from the smallest to the highest index). 
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If we add an extra initial vector u 3 = |en, Algorithm Q] terminates after 9 iterations with 
an invariant polytope of 48 • 2 vertices. The partial polytope norms are reported in Table 
0 The largest singular value is ay = 3.8190 • 10° and the smallest one is 09 = 1.5416 • 10 -5 , 
which demonstrate a relatively balanced shape. The computed Holder exponent is a^ 12 = 
12 - log 2 p(B 0 ) = 3.833483495658518 .... 


i 

4 5 6 7 

1 • lb Pi 

5.0579 1.5497 1.1597 1 


Table 5: The partial polytope norms computed by the modified algorithm for N = 12. 

6.2 The table of results for N < 20 

Proceeding this way we have computed the exact values of Holder exponents of Daubechies 
wavelets according to Table | 6 j 

We indicate the s.m.p., the extra initial vectors, the number of vertices of the final 
polytope (JfV), the number of iterations of Algorithm Q] (#its) and the Holder exponent a. 


N 

s.m.p. 

Extra vertices 

#its 

w 

a 

2 

B 0 

none 

1 

1 • 2 

0.55001... 

3 

Bo 

none 

3 

3-2 

1.08783... 

4 

Bo 

0 . 8 e 3 

2 

4-2 

1.61792... 

5 

B 0 and Bi 

0 . 1 e 4 

4 

8-2 

1.96896... 

6 

B 0 and Bi 

0 . 1 e 5 

5 

11 • 2 

2.18913... 

7 

B 0 and B\ 

0 . 1 e 5 

5 

12 • 2 

2.46040... 

8 

Bo and B\ 

0 . 1 e 7 

5 

18 • 2 

2.76081... 

9 

B 0 and B\ 

0.5e 8 

6 

24-2 

3.07361... 

10 

Bim 

0.5e 9 

10 

90-2 

3.36139... 

11 

B 0 and Bi 

0.5ei 0 

11 

75-2 

3.60346... 

12 

B 0 and B\ 

0.5en 

7 

48 • 2 

3.83348... 

13 

B 0 and B\ 

Cl 2 

18 

73-2 

4.07347... 

14 

B 0 and B\ 

0.5ei 3 , 0.25ei2 

15 

73-2 

4.31676... 

15 

B 4 0 Bf 

10 ~ 3 {e k }lU 

14 

376-2 

4.55611... 

16 

B 2 oB\ 

™- 2 W}lln 

13 

372 • 2 

4.78643... 

17 

Bq and Bi 

10 - 3 {e k }H u 

11 

480 • 2 

5.02444... 

18 

B 0 and B\ 

10 - 3 {e fc ULi 2 

13 

409 • 2 

5.23915... 

19 

B 0 and B\ 

10 - 3 {e,U 8 =13 

19 

1395 • 2 

5.46529... 

20 

B 0 and Bi 

10 - 3 {e fc U 9 =13 

24 

2480 • 2 

5.69116... 


Table 6 : Computed Holder exponent of Daubechies wavelets. 
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Appendix 1. Proof of Theorem HJ 

Necessity. If the algorithm terminates within finite time, then the products {Ih }[ =1 are 
dominant and their leading eigenvalues are unique and simple. This is shown in the same 
way as in the proof of Theorem 4 of El for r = 1. To prove that a is admissible, we take 
arbitrary i and j and denote by z the vertex of the final polytope P with the largest scalar 
product (v*,z). Since (y*. fPjz) = ([fl*] k v*,z) = (v*,z), all the points {ft^} fceN also provide 
the largest scalar product with the vector v*. Hence, they are all on the boundary of P, i.e., 
they are not absorbed in the algorithm. Consequently, the algorithm can terminate within 
finite time only if z is the leading eigenvector of n^, i.e., z = aqiq. Thus, the maximal scalar 
product ( v*,z ) over all z G P is attained at a unique vertex z = aqiq, where it is equal to 
(vj, otjVj) = ctj. Since i ^ j, it follows that 

sup (v*-, z) < aq. 

z£aiPi }00 

Thus, atqij < aq, which proves the admissibility of a. 

Sufficiency. Denote by O the set of products n^, i = 1 ,... , r, and of its cyclic permuta¬ 
tions. Since this is a set of dominant products for A, their leading eigenvectors { v^ | k = 
1 ,..., rii , i — 1,..., r} are all different up to normalization, i.e., they are all non-collinear. 
Indeed, if, say, vf ' i = A i/p , A ^ 0, then, replacing the products n* and by the corre¬ 
sponding cyclic permutations, it may be assumed that k — l — 1. However, in this case 
= n^ 2 uf } = v\ 1} for any k \, fc 2 - Therefore, the spectral radius of every product 
of the form is at lest one. By the dominance assumption, this product is a power 

of some product nefi. Taking now &q, /C 2 large enough and applying Lemma |2] first to the 
words a = n, 6 = and then to the words a = n, b = H;, we conclude that both nj, n, ; must 
be cyclic permutations of ft, which is impossible. Thus, all the leading eigenvectors 
are non-collinear. Hence, there is e > 0 such that for every x G \ {0} the ball of radius 
e||x|| centered at x may contain leading eigenvectors of at most one matrix from D. 

If the polytope algorithm with the initial roots aiT-Li ,..., a r T-L r does not converge, 
then there is an element of some root, say, aqui = aquj 1 * G aiHi and an infinite se¬ 
quence {^4fe fc }fceN; which is not periodic with period fti, and such that every vector = 
Ab k _i ■ ■ ■ A^aiVi is not absorbed in the algorithm. This implies that there is a constant 
Co > 0 such that Huftl > Co for all k. On the other hand, ||tifc|| < M ||o!i iq||, hence the 
compactness argument yields the existence of a limit point m ^ 0 of this sequence. Thus, 
for some subsequence, we have Uj k —> u as k —> oo. Let S > 0 be a small number to specify. 
Passing to a subsequence, it may be assumed that \\uj n — Uj k \\ < 5 for all k,n. Denote 
Gk = Aj k+1 _i • • • Aj k . We have Uj k+1 = GkUj k , k G N. Invoking the triangle inequality, we 
obtain 

||GfcM-u|| < \\G k {u - u jk )\\ + \\G k u jk - u jk \\ + \\uj k — w|| < M5 + 5 + 5 = (M + 2)6. 

Hence, Lemma [T] yields p(P k ) > 1 — C(d)M 1+ ^5 1 ^ d for all k. The dominance assumption 
implies that if 6 > 0 is small enough, then all G k must be powers of matrices from O. Each 


matrix from fl has a simple unique leading eigenvalue 1. Therefore, there is a function p(t) 
such that /i(t) —y 0 as t —> 0, and for every matrix Q which is a power of a matrix from 
Q the inequality || Qu — u || < t implies that there exists a leading eigenvector w of Q such 
that ||w — u|| < p{t). Thus, for every k G N there exists a leading eigenvector Wk of Gk 
such that ||Wk — w|| < /i((M + 2)5). For sufficiently small 5 we have /r((M + 2)5) < e||u||. 
Hence, for every k, the vector belongs to the ball of radius e||w|| centered at u. However, 
this ball may contain a leading eigenvector of at most one matrix from fi, say n. Therefore, 
all Gk, k G N, are powers of II and u is the leading eigenvector of IT. Thus, Uj k = tl Vk a.iVi 
for some pk G N. Clearly, n is a cyclic permutation of some Hp If j ^ 1, then assuming 
that n = Uj (the general case is considered in the same way), we have Uj k —> ai(v* ,v\)vj 
as k —> oo. Since the balancing vector a is admissible and (v*,vi) < q\ v it follows that 
a< atj. Therefore, the limit point u = aq(u*, V\)Vj = A ctjVj for some A G (0,1), 
is interior for the initial polytope co {aH). This means that for large k, the point Uj k will 
be absorbed in the algorithm, which contradicts to the assumption. Consider the last case, 
when n is a cyclic permutation of Hi. Since u is the leading eigenvector of n we see that 
u = /3v[ s ^ for some s = l,... ,ni and f3 G M. We assume f3 > 0, the case of negative (3 
is considered in the same way. If f3 < 1, then we again conclude that Uj k are absorbed 
in the algorithm for large k. If (3 > 1, then for the product n 0 = A dn • • • A ds , we have 
U 0 u = f3v i, and hence HoWf! —» (3v i. Therefore, ||/I _ 1 non A: ui — ui|| —$■ 0 as k —> oo. 
By Lemma [TJ this means that the spectral radius of the product HoH^ tends to (3 > 1 as 
k —y oo. The dominance assumption implies now that for every sufficiently large k, this 
product is a power of some II a G fl. Applying Lemma [2] to the words a = Il„, b = IT, we see 
that n a is a cyclic permutation of n. In particular, |II a | = |n| = n. Therefore, the length 
|n 0 W| = (n — s + 1) + kn = {k + l)n — (s — 1) is divisible by |II a | = n. Hence, the number 
(s — 1) is divisible by n, which is impossible, because s < n. This completes the proof. □ 

Appendix 2: the Butterfly matrices. 

We report here the matrices of the Butterfly scheme for oo = A, First the 17 x 17 family A, 
then the 11x11 family B and finally the 6x6 family C. 
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The four matrices Ai, A 2 , A%, A 4 : 


34169 

4104 

_ 1319 
2052 

5081 

1026 

2515 

1026 

_ 11509 
5130 

_ 42173 
10260 

163009 

20520 


33437 

5130 

132199 

20520 


3617 

2565 


21365 

4104 


773 
' 228 


1283 

456 


181 

1140 

679 

380 

1411 

1140 

1021 

1140 


_ 523 
855 

3991 

2280 

4129 

1710 


9313 

6840 


_ 4225 
” 1368 

4133 

1368 


4643 

1368 

4577 

1368 

112 

171 

365 

342 

167 

180 

784 

285 

5249 

1140 


509 

1140 

7453 

1140 


193 

1368 


3062 

513 

475 

108 


2677 

4104 

4501 

2052 

3661 

1026 

3935 

1026 

4391 

5130 

21967 

10260 

66761 

20520 


3284 

2565 

35951 

20520 


151 

270 


403 

216 


7531 

4104 

5893 

4104 


10537 

4104 


991 

5130 

7168 

2565 

13723 

10260 

533 

270 

2971 

2565 


2863 

2052 

205 

513 

5068 

2565 


7219 

4104 


120 

19 


5017 

3420 

1061 

1710 


688 

171 

11701 

6840 

3253 

3420 


4751 
' 6840 


7691 

684 


16343 

1368 


4319 

684 


2377 

684 

_ 824 
285 


1829 

1710 

1931 

380 


125431 

6840 


9157 

1368 


16343 

1368 

1571 

342 

_ 2345 
1368 


4319 

684 


824 
' 285 


1931 

380 

1829 

1710 

22037 

2280 


2081 

570 


11009 

4104 

2813 

2052 

1132 

513 

11663 

4104 

6641 

2052 


20821 

5130 


17159 

5130 

13063 

5130 

5057 

5130 


2605 
' 2052 

734 

513 


19553 

2052 


4561 
' 2052 


9409 

4104 

1183 

1026 

20803 

4104 

2549 

513 

8057 

2052 


11303 

10260 

2573 

2565 

7679 

2565 


_ 8411 
5130 

30991 

10260 


121739 

20520 


10501 

4104 

4561 

2052 


20803 

4104 

1183 

1026 

9409 

4104 

11303 

10260 


8057 

2052 

7679 

2565 

2573 

2565 


30991 

10260 

_ 8411 
5130 

66763 

20520 


3703 

1026 


2813 

2052 


6641 

2052 

11663 

4104 


20821 

5130 

5057 

5130 


734 

513 

2605 

2052 

3841 

20520 


31475 

4104 


3673 

1368 


2281 

456 

449 

342 


14287 

3420 

6337 

1710 

1333 

3420 


674 

855 

139 

684 

4843 

1368 


619 

1710 


2399 8869 

684 1368 

3851 785 

1368 456 

1643 _ 4919 

1368 1368 

1573 3469 

684 456 

409 6971 

684 1368 

1525 677 

1368 152 

1309 1657 

684 1140 

1189 599 

855 228 

3181 401 

1710 285 

107 _ 10439 

228 3420 

14827 2009 

6840 380 

29 1189 

57 684 

3233 47 

3420 380 

_ 179 6845 

760 684 

2231 109 

1710 380 

1859 29269 

570 6840 

6211 13183 

3420 6840 


397 286 

1368 57 

397 3673 

1368 1368 

1603 1033 

1368 456 

2275 449 

1368 342 

2275 2281 

1368 456 

1603 1939 

1368 342 

2771 14287 

3420 3420 

673 4241 

855 1710 

2771 1097 

3420 342 

1663 1333 

1140 3420 

1663 6337 

1140 1710 

8 4843 

19 1368 

2713 139 

1140 684 

2713 674 

1140 855 

8 63959 

19 6840 

23827 619 

6840 1710 

23827 13529 

6840 2280 


27103 18013 

4104 4104 

15553 7957 

4104 1026 

15629 7639 

4104 2052 

30277 1795 

4104 4104 

14153 14935 

4104 4104 

19439 _ 391 

4104 513 

61 7141 

540 2565 

2255 3863 

2052 10260 

17719 7006 

10260 2565 

_ 25 _ 12499 

1026 10260 

14051 118553 

2052 20520 

20527 6691 

10260 2052 

14779 9307 

10260 5130 

65351 _ 40487 

10260 20520 

22733 1093 

10260 2052 

77327 _ 52403 

20520 20520 

22907 84623 

20520 20520 


10133 6077 

2052 4104 

_ 6077 10133 

4104 2052 

5915 5657 

4104 1026 

2411 20279 

2052 4104 

20279 2411 

4104 2052 

5657 5915 

1026 4104 

21263 953 

10260 2052 

_ 1849 _ 1849 

10260 10260 

953 21263 

2052 10260 

16613 8501 

5130 2052 

8501 16613 

2052 5130 

_ 68621 24209 

20520 4104 

6062 40547 

2565 10260 

40547 6062 

10260 2565 

24209 68621 

4104 20520 

4469 193087 

10260 20520 

193087 4469 

20520 10260 


7727 _ 1870 

4104 513 

9767 _ 8171 

4104 4104 

5251 9871 

4104 4104 

1897 _ 2441 

4104 1026 

7043 1547 

4104 2052 

5069 4825 

4104 4104 

1139 1963 

2565 10260 

965 13877 

1026 10260 

5029 8209 

10260 10260 

2981 7787 

1026 2565 

2269 9701 

2052 4104 

26377 9703 

10260 10260 

6526 _ 4882 

2565 2565 

5029 15133 

10260 4104 

13849 19073 

5130 10260 

48463 32657 

20520 10260 

4067 12139 

20520 10260 


2209 13 

1026 216 

_ 8237 319 

4104 4104 

20095 13297 

4104 4104 

101 7321 

27 4104 

5789 5665 

4104 4104 

1601 11569 

2052 4104 

6568 4877 

2565 2052 

_ 5416 _ 24079 

2565 10260 

1231 7367 

2052 10260 

61 _ 1081 
135 10260 

2041 2455 

2052 2052 

78761 2923 

20520 2565 

1303 35 

2565 54 

3188 1259 

2565 10260 

139 8251 

4104 2565 

17584 79721 

2565 20520 

106477 9343 

20520 20520 


_ 4225 _ 6887 

2052 4104 

_ 449 _ 17 

216 54 

15877 7631 

4104 2052 

_ 299 _ 7655 

1026 4104 

4181 14425 

2052 4104 

12691 4303 

4104 1026 

8447 _ 977 

10260 5130 

1199 3193 

513 2565 

_ 3826 _ 20899 

2565 10260 

20039 2999 

10260 2565 

39097 32347 

20520 20520 

13943 _ 793 

5130 1026 

_ 2609 _ 2063 

2052 5130 

_ 2051 _ 4313 

1080 1080 

415 160 

2052 513 

8377 19687 

5130 20520 

_ 7481 _ 54053 

5130 20520 


319 _ 8237 

4104 4104 

13 2209 

216 1026 

11569 1601 

4104 2052 

5665 _ 5789 

4104 4104 

7321 101 

4104 27 

13297 20095 

4104 4104 

7367 1231 

10260 2052 

24079 _ 5416 

10260 2565 

_ 4877 _ 6568 

2052 2565 

2455 2041 

2052 2052 

1081 61 
10260 135 

_ 8251 _ 139 

2565 4104 

1259 3188 

10260 2565 

35 1303 

54 2565 

_ 2923 78761 

2565 20520 

9343 106477 

20520 20520 

79721 17584 

20520 2565 


1099 887 \ 

456 1368 

47 3911 

171 1368 

_ 443 _ 2347 

228 1368 

5005 _ 2347 

1368 1368 

1429 3911 

456 1368 

523 887 

171 1368 

3871 _ 191 

3420 855 

3107 _ 191 

3420 855 

53 _ 106 

171 171 

_ 115 _ 2429 

342 3420 

2169 _ 1658 

760 855 

5467 775 

3420 342 

761 775 

1140 342 

6557 _ 1658 

1368 855 

1193 _ 2429 

1140 3420 

69 637 

760 6840 

7027 _ 7129 , 

6840 6840 / 


5005 443 \ 

1368 228 

_ 443 5005 

228 1368 

47 _ 1429 

171 456 

1099 523 

456 171 

523 1099 

171 456 

_ 1429 47 

456 171 

3871 _ 3107 

3420 3420 

53 53 

171 171 

_ 3107 3871 

3420 3420 

_ 761 _ 5467 

1140 3420 

_ 5467 _ 761 

3420 1140 

_ 2169 6557 

760 1368 

_ 115 1193 

342 1140 

1193 _ 115 

1140 342 

6557 _ 2169 

1368 760 

21281 1487 

6840 285 

1487 _ 21281 , 

285 6840 / 





















13997 15049 

1026 4104 

22307 54187 

4104 4104 

_ 1319 _ 41957 

2052 4104 

34169 _ 43913 

4104 4104 

17155 65953 

4104 4104 

16895 29369 

2052 4104 

11509 1003 

5130 10260 

2515 _ 929 

1026 270 

_ 5081 889 

1026 513 

163009 _ 56513 

20520 5130 

42173 8518 

10260 2565 

_ 3617 11347 

2565 5130 

132199 92413 

20520 5130 

33437 23953 

5130 10260 

51557 _ 51289 

5130 10260 

21365 44899 

4104 20520 

79909 30997 

20520 4104 


7531 3062 

4104 513 

8251 _ 15943 

4104 4104 

14623 4501 

4104 2052 

10537 _ 2677 

4104 4104 

133 1415 

216 4104 

5893 475 

4104 108 

13723 4391 

10260 5130 

7168 3935 

2565 1026 

991 3661 

5130 1026 

2971 66761 

2565 20520 

_ 533 21967 

270 10260 

_ 5068 _ 151 

2565 270 

205 35951 

513 20520 

2863 _ 3284 

2052 2565 

3443 _ 12344 

10260 2565 

_ 7219 403 

4104 216 

8299 78761 

4104 20520 


_ 4225 89 

1368 342 

10201 131 

1368 456 

4577 863 

1368 171 

4643 _ 1283 

1368 456 

_ 9643 1475 

1368 1368 

4133 773 

1368 228 

167 1411 

180 1140 

365 679 

342 380 

112 181 
171 1140 

5249 _ 11933 

1140 6840 

784 1021 

285 1140 

43 4129 

285 1710 

_ 7453 3991 

1140 2280 

509 523 

1140 855 

75 21 

76 95 

193 9313 

1368 6840 

42259 2113 

6840 1368 


785 3851 

456 1368 

8869 2399 

1368 684 

677 1525 

152 1368 

6971 409 

1368 684 

3469 1573 

456 684 

4919 1643 

1368 1368 

401 3181 

285 1710 

599 1189 

228 855 

1657 1309 

1140 684 

2009 14827 

380 6840 

10439 107 

3420 228 

109 2231 

380 1710 

6845 _ 179 

684 760 

47 3233 

380 3420 

1189 29 

684 57 

13183 6211 

6840 3420 

29269 1859 

6840 570 


1453 11939 

152 1368 

11939 1453 

1368 152 

743 5315 

76 1368 

9425 5 _ 

684 152 

_ 5 _ 9425 

152 684 

5315 743 

1368 76 

59581 _ 5923 

3420 684 

6187 6187 

684 684 

5923 59581 

684 3420 

37273 _ 1471 

6840 380 

1471 37273 

380 6840 

30023 12133 

6840 2280 

6499 433 

2280 3420 

433 6499 

3420 2280 

12133 30023 

2280 6840 

_ 99 32563 

38 6840 

32563 _ 99 

6840 38 


307 1985 

152 1368 

2507 439 

1368 152 

_ 21 2693 

8 684 

1921 _ 431 

1368 “ 76 " 

27 _ 95 

76 72 

1091 385 

684 152 

2881 21659 

3420 3420 

5411 2201 

3420 684 

1561 1937 

3420 684 

9607 11699 

3420 2280 

6017 56 

2280 45 

24793 4019 

6840 2280 

145 21673 

228 6840 

161 187 

1368 95 

_ 187 23801 

40 6840 

_ 511 _ 200 

152 171 

59 _ 353 

18 760 


439 2507 

152 1368 

1985 307 

1368 152 

385 1091 

152 684 

_ 95 27 

72 76 

_ 431 1921 

76 1368 

2693 _ 21 

684 8 

1937 1561 

684 3420 

2201 _ 5411 

684 3420 

21659 2881 

3420 3420 

56 6017 

45 2280 

11699 9607 

2280 3420 

23801 _ 187 

6840 40 

187 161 

95 1368 

21673 145 

6840 228 

4019 24793 

2280 6840 

_ 353 59 

760 18 

_ 200 _ 511 

171 152 


6071 73 

1368 342 

5017 73 

1368 342 

4721 859 

1368 1368 

545 3317 

1368 1368 

227 3317 

36 1368 

3319 859 

684 1368 

_ 5729 _ 143 

1710 1710 

_ 8267 479 

1710 342 

3221 _ 143 

342 1710 

889 2113 

342 1368 

1237 _ 2113 

360 1368 

26437 _ 1183 

6840 3420 

851 17263 

855 6840 

865 17263 

1368 6840 

25979 _ 1183 

6840 3420 

16873 4117 

6840 6840 

1997 4117 

3420 6840 


3299 

1368 


7957 

1026 

18013 

4104 

391 

513 


15553 

4104 

27103 

4104 


6887 

4104 

4303 

1026 


_ 449 
216 

_ 4225 
2052 

12691 

4104 


8171 

4104 

1870 

513 


9767 

4104 

7727 

4104 


3911 

1368 

887 

1368 

887 

1368 

3911 

1368 


47 

171 

1099 

456 


1429 

456 


980 

171 

1061 

1710 

5017 

3420 


7639 

2052 

7006 

2565 

3863 

10260 


15629 

4104 

17719 

10260 

2255 

2052 


7631 

2052 

20899 

10260 

3193 

2565 


15877 

4104 


1199 

513 


9871 

4104 


5251 

4104 


2347 

1368 


191 

855 


443 

228 


3107 

3420 


1481 

360 


3253 

3420 


118553 

20520 

_ 12499 
10260 

1093 

2052 


14051 

2052 


22733 

10260 


32347 

20520 

2999 

2565 

160 

513 


39097 

20520 

20039 

10260 


9701 

4104 


2269 

2052 


13849 

5130 


1658 

855 

2429 
' 3420 

2429 

3420 


2169 

760 


18733 

3420 


9307 

5130 

6691 

2052 


14779 
‘10260 

20527 

10260 


2063 
' 5130 

793 

1026 


13943 

5130 


775 

342 


775 

342 


5467 

3420 


2611 
' 2280 


52403 

20520 


77327 
' 20520 


19687 

20520 


8377 

5130 


32657 

10260 


48463 

20520 


637 

6840 


69 

760 


5017 6221 2933 3911 989 517 2003 463 385 \ 

1368 1368 1368 1368 684 684 1368 152 152 


5017 

1368 

6071 
" 1368 


6221 

1368 

2933 
" 1368 


545 

1368 

4721 
‘ 1368 


5729 

1710 

1237 

360 

889 

342 


865 

1368 

851 

855 


16873 

6840 


1289 

171 

3209 

1368 

2705 
" 1368 


1031 

380 

1931 

1368 

4187 
" 3420 


15277 

6840 

3169 

3420 

9743 

6840 


6241 

6840 


2933 

1368 

6221 

1368 

2705 

1368 

3209 

1368 

1289 

171 


1031 

380 


4187 

3420 

1931 

1368 


3169 

3420 


1507 

1368 


7913 

3420 


3911 

1368 

2003 
' 1368 


2431 
' 1368 

6145 

1368 


3913 

1140 


1289 

1710 


667 

6840 


4291 

6840 


517 

684 

989 

684 

4613 

1368 

1301 

1368 

2113 
' 1368 

3997 

1368 

241 

228 


3797 

1140 

7879 

6840 


569 

6840 

683 

1368 


23003 

6840 


2003 

1368 

3911 

1368 

6145 

1368 

2431 

1368 


3913 

1140 


1289 

1710 


43 

1710 

667 

6840 

1307 

1368 


6919 

1710 


463 

152 

385 

152 

385 

152 

463 

152 

39 

76 

39 


677 

285 

26 

57 

2569 

2280 


1529 

2280 

1529 

2280 


1091 

570 


463 

152 

39 

76 

39 

76 

463 

152 

385 

152 

677 
' 285 


1354 

285 

2569 

2280 


1529 

2280 


2569 

2280 


2803 

2280 









the four matrices Bi, B 2 , B 3 ,- 84 : 


/ - 


b 2 = 


b 3 = 


s 4 = 


245695 

30276 

27385 

7569 

2425 

841 


825575 

30276 


6728 

738031 

20184 

1401 

29 

4781 

174 

1650899 

121104 


121104 

467615 

40368 

148535 

40368 

1148275 


779423 

26912 


285081 

26912 


3977 

232 

13555 

696 

79195 

15138 


14795 

5046 

10687 

1682 


37685 

15138 

49977 

6728 


2415 

116 

6433 

348 


1024085 

121104 


121104 
200725 


36485 

40368 


121104 

616655 


47955 

26912 


_ 3977 
232 

13555 

696 


159509 

121104 


53651 
'13456 


403680 

1481557 

403680 

21023 

3480 

12787 

3480 


68083 

60552 

40939 

60552 

49109 

60552 


659713 

60552 


201840 

1203897 

67280 

32659 

1740 

19261 

1740 


121104 

383841 

13456 

55715 

40368 

1918919 

121104 

2091859 

121104 

4879959 

134560 

1103137 

134560 

27188873 


178027 

3480 


940283 
' 121104 


121104 

1294319 


103615 

13456 


121104 

3060847 

134560 

1160719 

134560 

21731537 


178027 

3480 


121104 

231455 

121104 

462383 

40368 

113975 

40368 

908875 

121104 

686915 

121104 

366567 

26912 

486971 

80736 

1852889 


13987 

696 


499045 
‘ 121104 


17775 

13456 

547705 

121104 


_ 465245 
121104 

263909 

26912 

253961 

80736 


8521 

696 


_ 46215 
13456 

575405 

121104 

183455 

121104 

252649 

121104 

245095 

121104 

120045 

26912 

424769 

80736 

342619 

80736 

1055 

696 

893 

232 


293725 

121104 

169985 

121104 

344065 

40368 

54145 

13456 

99575 

121104 

679265 

121104 

130045 

26912 

454103 

80736 

964451 


32845 

30276 

6745 

30276 

22085 
"30276 

950 

7569 

9925 

7569 


2486 

2523 

163 

1682 


152425 

40368 

150395 

40368 

967175 

121104 

118385 

121104 

144065 

121104 


80736 

121883 

80736 

4775 

696 

2075 

232 


107765 

121104 


1945 

13456 

51305 

40368 

430865 

121104 

31715 

121104 

153665 

26912 

124463 

26912 

556927 


1985 

696 


8521 

696 


121104 

325415 

121104 

154785 

13456 

57715 

40368 

740095 

121104 

692831 

121104 

314715 

26912 

258103 

80736 


13987 

696 


465355 

121104 

384635 

121104 

165345 

13456 

184385 

40368 

149515 

121104 

1055135 

121104 

185503 

26912 

461357 

80736 

1588793 


13207 

696 


679825 

30276 


84875 

5046 

29950 

2523 


24497 

3364 

_ 81551 
3364 

154865 

10092 

1401 

29 

4781 

174 

71365 
' 15138 


57415 

5046 

24185 

3364 

296593 

15138 

254315 

15138 

4762 

841 


2415 

116 

6433 

348 


8855 

7569 


43595 

10092 

7235 

10092 

9085 

15138 


10645 

1682 


192245 

40368 

192535 

40368 

_ 251765 
121104 

203869 

121104 

134675 

121104 

78545 

121104 

300223 

26912 

287819 

80736 

918185 

80736 

1055 

696 

893 

232 

10381 

30276 


25265 

10092 

31055 

10092 


39545 

30276 

15721 

3364 


1795 

348 


555283 

60552 


60552 

381797 


248029 

60552 


8593 

67280 

1697971 

201840 

353907 

67280 

32659 

1740 

19261 

1740 

13805 

30276 


44915 

10092 

8285 

1682 

262085 

30276 

170995 


32851 

6728 


1795 

348 


121104 

553669 

121104 

94085 

121104 

139145 


303419 

134560 

734231 

403680 

236413 

403680 

21023 

3480 

12787 

3480 

38465 

30276 


12851 

5046 

1525 

1682 


29733 

6728 


338935 

30276 


46190 

7569 


15529 

6728 

247579 

20184 

138401 

20184 

2039 


441205 
' 121104 


66955 

40368 

16795 

40368 

327335 

121104 

186655 

121104 

62849 

26912 

30097 

26912 

92615 

80736 

1145 

232 

1985 

696 


33280 

2523 


20184 

361817 


51605 

15138 

25345 

30276 

2165 

7569 


30625 

2523 


20184 

326087 


7421 

348 


12635 

7569 

5675 

7569 

22400 

7569 

19955 

30276 

2305 

5046 


2953 

2523 

4297 

10092 

205 

87 

1313 

348 


5 

696 


1075 

696 

1075 

696 

235 

696 

287 

464 

1085 

464 


842065 

121104 


695 

13456 


375917 

26912 

219041 

80736 


13207 

696 


220075 

40368 


69905 

40368 

94545 

13456 

81845 

13456 


436517 

80736 


4417 

232 


18395 
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and the four matrices C’i, C 2 , C :i , C 4 : 
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c 2 


c 3 


c 4 
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0 10 0 J 0 
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12 16 12 


\ 0 0 

( ° 3 

1 3 
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0 
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0 

3 

4 
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4 

0 
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12 

25 

12 
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16 
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1 
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Ao _ Zo q 
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1 J 
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